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Abstract 

A feasible model is introduced that manifests phenomena intrinsic to iterative complex 
analytic maps (such as the Mandelbrot set and Julia sets) . The system is composed of two 
coupled alternately excited oscillators (or self-sustained oscillators). The idea is based on 
a turn- by-turn transfer of the excitation from one subsystem to another (S.P. Kuznetsov, 
Phys. Rev. Lett. 95 , 2005, 144101) accompanied with appropriate nonlinear transformation 
of the complex amplitude of the oscillations in the course of the process. Analytic and 
numerical studies are performed. Special attention is paid to an analysis of the violation of 
the applicability of the slow amplitude method with the decrease in the ratio of the period of 
the excitation transfer to the basic period of the oscillations. The main effect is the rotation 
of the Mandelbrot-like set in the complex parameter plane; one more effect is the destruction 
of subtle small-scale fractal structure of the set due to the presence of non-analytic terms in 
the complex amplitude equations. 

PACS.Q5.45.-& 
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1 Introduction 



One special chapter of nonlinear dynamics elaborated extensively by mathematicians consists in 
a study of iterative maps denned by analytic functions of complex variable. A classic object is a 
complex quadratic map pQ 

z n+1 = c + z 2 n . (1) 

At c = the behavior of the iterations is rather evident: for \zq\ > 1 the result diverges to 
infinity, and for \z Q \ < 1 one observes residence of the variable z in a bounded part of the complex 
plane, and indeed convergence to 0. The border between these two kinds of behavior is the unit 
circle \zq\ = 1. At other values of the complex parameter c a set in the complex plane z, the 
Julia set, separating the bounded from divergent behavior appears to be rather complicated and 
nontrivial (fractal). See examples in Fig. [TJ 
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Figure 1: Mandelbrot set (a) and Julia sets for the complex quadratic map ([I]) at several values of 
the parameter c = 0.0 (b), c = 0.2-0.3? (c), c = -0.4 (d), c = -0.8 (e), c = -0.1 + 0.75z (f). Gray 
areas correspond to periodic dynamics in a bounded domain of the complex variable z (the periods 
are marked by figures). Black designates bounded chaotic dynamics, and white corresponds to 
escape of the iterations to infinity. 
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Alternatively, we can fix the initial condition zq = and consider how the behavior of the 
iterations depends on the complex parameter c. Then, at some values of c the iterations escape to 
infinity, and at others they stay in a bounded domain. The set on the complex plane c associated 
with the latter situation is called the Mandelbrot set. See Fig. [Ha). Note that the bounded 
dynamics may be periodic, quasiperiodic or chaotic. The domain of periodic behavior resembles 
a cactus and consists of a set of roundish formations touching each other and placed around one 
large area having the form of a cardioid. Figures 1, 2, 3 on the plot designate periods of dynamics 
observed in several basic leaves of the "cactus" . The dots in Fig. UJa) indicate parameter values, 
for which the Julia sets are shown in panels (b)-(f). 

A challenging problem is the realization of the phenomena intrinsic to the complex analytic 
maps in physical systems (see, e.g. [2]). One successful attempt to implement this dynamics relates 
to a symmetric system of coupled maps for two real variables or to a system of two periodically 
driven coupled nonlinear oscillators [3J. In this context an interpretation of the Mandelbrot set 
was suggested, as a domain of generalized partial synchronization [I], and some aspects of the 
feasibility of the phenomena of complex dynamics in autonomous flow systems were studied [5]. 
An electronic analog device was designed simulating the dynamics of two coupled quadratic maps, 
in which the Mandelbrot-like set was observed for the first time in a physical experiment [6j. 

In the present paper, we suggest an alternative approach that gives an opportunity to organize 
dynamics similar to that in the complex analytic maps. The main idea is based on interpretation 
of a complex amplitude of an oscillatory process as a complex variable. Two alternately excited 
oscillators are used to pass the excitation between each other with transformation of the complex 
amplitude corresponding approximately to the complex quadratic map. Earlier a similar idea was 
applied to construct physical systems delivering realistic examples for some well-known abstract 
concepts and phenomena of the nonlinear science, including the Bernoulli map, a Smale- Williams 
attractor, Arnold's cat map, and a robust strange nonchaotic attractor [TJ, El El Ell ITT] - 

In Sec. 2 we introduce a system of coupled oscillators alternately activated by means of mod- 
ulation of the dissipation parameter accompanied with a turn-by-turn transfer of the excitation 
between the subsystems. In Sec. 3 we derive the shortened equations for the system using the 
complex amplitude approach. Then, we undertake an approximate analytic solution of the equa- 
tions. The result is the complex quadratic mapping, which represents the Poincare map for the 
system and governs the state evolution on one period of the parameter modulation. We present 
and compare pictures on the parameter plane and the phase space portraits for the amplitude 
equations and those for the approximate analytic map (the Mandelbrot set and the Julia sets). In 
Sec. 4 we turn again to the original coupled oscillator equations and consider results of numerical 
studies to observe and discuss deviations from the slow-amplitude approach, which become rele- 
vant with decrease of a parameter N representing a ratio of the modulation period to the basic 
period of oscillations. We find that the most notable effect consists of rotation of the Mandelbrot- 
like set in the complex parameter plane while its visible fractal-like structure persists. In addition, 
we reveal the gradual destruction of subtle details of this fractal structure (starting from smaller 
scales) as the parameter N departs from the range of applicability of the method of slow ampli- 
tudes. We relate this phenomenon to non-resonance complex conjugate terms in the equations of 
the oscillators (rewritten in terms of the complex variables), which are negligible in the large N 
asymptotics. In Sec. 5 we modify original system in such way, that it becomes a system of coupled 
self- sustained oscillators. We show that in the parameter plane of this system the Mandelbrot-like 
set can preserve. Such system of two coupled elements, based on van der Pol oscillators, is more 
convenient for electronic device construction. Moreover, in special case it is equivqlent to the 
system with a Smale- Williams attractor [ID] . Thus, possibility of coexistence of the phenomena 
of complex analytic dynamics and the phenomenon of hyperbolic chaos is manifested. 
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2 The basic model system 



In the theory of oscillations and waves, the method of slow amplitudes is well-known. An os- 
cillatory process possessing a basic frequency ujq is attributed with the slow complex ampli- 
tude A(t) as follows: x(t) = A(t)e iuJot + A*(t)e~ iuJot = 2Re[A(t)e iuK)t ], where asterisk designates 
complex conjugate. Transformation of this signal by means of a quadratic nonlinearity yields 
y{t) = [x(t)} 2 = 2\A(t)\ 2 + 2Re{[A(t)] 2 e 2tuJot } . As we may see, the doubled frequency component 
has the complex amplitude equal to the squared complex amplitude of the original signal. Such a 
transformation will be a main element of the approach we develop in the realization of phenomena 
of complex analytic dynamics. 

Let us consider a system of two coupled non-autonomous oscillators 



where x and y represent generalized coordinates for the first and the second oscillator, respectively, 
and F, 7, A, ip, e are parameters. The first oscillator has a characteristic frequency cuq, and the 
second one has a frequency twice as large. Parameters controlling dissipation in both oscillators 
vary slowly in counter-phase. The period of modulation T = 2ir/Q is assumed equal to an integer 
number N of periods of the basic oscillations 2-k/loq. 

One can describe the dynamics in terms of discrete time by means of a Poincare map. Given a 
vector x n = {x(t n ),u(t n ),y(t n ), v(t n )} = {x(t n ), x(t n )/u , y(t n ), y(t n )/2u } as the state of the sys- 
tem at t n = nT + to, from the solution of the differential equations ([2]) with the initial condition x n , 
we get a new vector x n+1 at t n+1 = (n + l)T +t . As it is determined uniquely by x n , we may intro- 
duce a function that maps the four-dimensional space {x,u,y,v} into itself: x n+1 = T(x n ). This 
Poincare map appears due to evolution determined by the differential equations with smooth and 
bounded right-hand parts in a finite domain of variables {x, u, y, v}. In accordance with theorems 
of existence, uniqueness, continuity, and differentiability of solutions of differential equations, the 
map T is a diffeomorphism, a one-to-one differentiable map of class C°°. 

Let us consider qualitatively how the system ([2]) behaves. The constant 7 is assumed to be 
positive, less than 1. Being on average positive, on a certain part of the modulation period, 
however, the dissipation parameter F{^/ + sin fit) for each oscillator becomes negative. On that 
interval, the oscillator is active (the oscillations grow). On the remaining part of the modulation 
period, it is passive (the oscillations decay). Let us assume that at the beginning of the active 
stage of the second oscillator, the first one oscillates with complex amplitude A: x(t) cx ~Re(Ae lulot ) . 
Then, the "germ" for excitation of the second oscillator (see the right-hand part of the second 
equation) is the second harmonic component produced by the nonlinear quadratic transformation 
of the signal from the first oscillator, Ke(A 2 e 2luJot ). The complex amplitude of the second oscillator 
on its active stage will be proportional to A 2 . Mixing with the auxiliary signal (see the right-hand 
part of the first equation) gives rise to a frequency component u with amplitude proportional 
to A 2 . Then, the sum of this component with an additional oscillatory term of frequency Uq, 
amplitude A and phase (p, produces a "germ" for the excitation of the first oscillator. Its complex 
amplitude is proportional to A 2 plus a complex constant. Hence, the stroboscopic Poincare map 
reduced in a certain approximation to a two-dimensional map which, expressed in terms of complex 
amplitudes, corresponds to the complex quadratic map. The complex parameter is given in terms 
of the modulus A and the argument <p. Dependent on this parameter and initial conditions, it 
may happen that the solution for the equations of coupled non-autonomous oscillators is either 
bounded or escapes to infinity. A domain on the complex parameter plane Xe lLp where bounded 
attractors persist will be the analog of the Mandelbrot set. Sets in the phase space separating 
initial conditions of bounded and unbounded motions will be analogs of the Julia sets. 



x + uj 2 x + F ■ (7 + sinf2t)x = eysmuot + Xsm(coot + y?), 
y + (2uo) 2 y + F ■ (7 — sin Qt)y = ex 2 , 
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3 Amplitude equations, slow-amplitude asymptotics and 
derivation of an approximate Poincare map 

To derive equations for complex amplitudes in the system of two coupled oscillators we set 

x = Ae iuJot + A*e~ iw °\ y = Be 2ioJot + B*e~ 2iuJot (3) 

with additional conditions 

Ae iwot + A*e~ iuJot = 0, Be 2iwot + B*e' 2iuJot = 0. (4) 

Substitution into the equations (|2J) yields 

2iuj Ae iwot + ioo (Ae iu)ot - A*e~ iu;ot )F(^ + sin fit) = 

e(Be 2iuJot + B*e- 2iuJot ) sin u t + X sm(u t + (p), (5) 
4iu Be 2iuJ0t + 2iu ( y Be 2iuJ ° t - B*e~ 2iuot )F( 1 - sin fit) = e(Ae iuot + A*e~ iuJot ) 2 , 

or, after multiplication of the first and the second equations by {2iujQe luJot )~ l and (4icj e 2 ' lU}ot )~ 1 , 
respectively, 

A + (F/2)(7 + sin fit) (A - A*e~ 2iuJot ) = 

-(e/2u )(Be 2iuiat + B*e- 2iujQt ){l - e ~ 2iuiat ) - (A/2o; )(e^ - e" 2 ^'-^), (6) 
B + (F/2)(7 - sinfit)( J B - B*e~^ ot ) = (e/Muq)(A + A* e~ 2iuJot ) 2 . 

The relations (Q are equivalent to the original equations (j2J) completely, although written in terms 
of the complex amplitudes A and £>. 

We now assume that the variation of the amplitudes A and B in time is slow, and that one 
can neglect the fast oscillating terms in ($5$). Formally, they are excluded by means of averaging 
the equations over a period 2tiu)q 1 . This approximation is justified in the case of a large frequency 
ratio N = uq/Q ^ 1. In this way, we arrive at the shortened slow-amplitude equations 



A + (F/2)(7 + sinfit)A = eB/(Au ) - Ae^/(4w ), 
B + (F/2)( 7 - sinfii)^ = eA 2 /(4iu ). 



■2,,.:., (7) 



With some additional assumptions, it is possible to construct analytically an approximate Poincare 
map for these equations. 

Let us start with consideration of the subsidiary differential equation 

w + g'(t)w = K(t)e m , (8) 

where g and / are some smooth real-valued functions. We represent a solution as 

w(t) = C(t)e- 9{t \ (9) 

where the time-dependent coefficient C(t) satisfies 

C = K(t)e f ® +a ®. (10) 

Consider a time interval containing a single maximum of the function f(t) + g(t) at t = to'. 

f (to) + g'(t ) = 0, /" (to) + /(to) = -P < o. (11) 

Integral of the equation (flQl) over the time interval can be evaluated with a help of the Laplace 
method. It yields 



K{t)e f{t)+9{t) dt « K(t )e f(to)+9(to) / e~ m 12 dt = K(t )e f{to)+ ^ to) y f 2^. (12) 



Let the solution (Q be zero to the left side from the maximum, then 



°® ~ { if (t ) v^7^e^Hf > t + h, (13) 

and 

/ .x _ / 0,t<t o -h, 

W[ ) ~ \ K(t )^/Pe^ + ^e-^\ t>t + h. 1 } 

Here /i may be thought as a characteristic width of the "hump" of the function / + g. 

Now we return to the amplitude equations (JTj). Let us consider separately two parts of one 
period of modulation. 

Firstly we consider the part containing the time interval of activity of the first oscillator. Here 
its amplitude is relatively large, while the second oscillator is passive and possesses very small 
amplitude. Here we may include only the effect of the first oscillator on the second one and 
neglect the backward coupling, i.e. regard the coupling as unidirectional. 

Excluding the right-hand part of the first equation, we write down for the slow amplitude of 
the first oscillator 

A(t) = C A exp 



p 

— (jt - fT 1 COS fit) 



(15) 



Let us define the value t n as associated with a maximum of \A\. It satisfies the relations 



sinfit n = —7, cosfit n = y 1 — 7 2 > 0. (16) 
Solving (FIBl) in respect to t n , substitute it into (IT51) and designate A n = A(t n ). Then, 



C A = A n exp 



il harctan-^ + v/W 



(17) 



For the second oscillator we write down 



B + -(7 - smQt)B = -^C 2 A exp [-F(jt - fi" 1 cos fit)] . (18) 

To solve this equation we use the result (fUj) obtained for the subsidiary problem (jSJ). To do this, 
we have to set 

0a-b(*) = (i72)( 7 -sinfit), /a-b(*) = -F( 7 t - fi- 1 cosfit), 

K A ^ B = [e/(Mu )}C 2 A liyj 

and 

9A^ B (t) + f A -> B {t) = -(l/2)F 7 t + (3/2)^-^08 fit, 

0A-b(*) + /a-b(*) = -(1/2)^7 - (3/2)Fsinfit, (20) 
9A^ B {t) + fl+B® = -(3/2)Fficosfit. 

From the condition of vanishing first derivative and the requirement for the second derivative to 
be negative, we determine the time instant t A ^B, the neighborhood of which corresponds to the 
excitation of the second oscillator by the first one. It satisfies 



sinfit^^s = —7/3, cosfit^^s = y 1 — 7 2 /9. (21) 

Here 

g' A -> B (tA^ B ) + /a->bIV*b) = -Pa^e = -(3/2)Fficosfit^ B = -(3/2)Ffiv/l-7 2 /9. (22) 
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Applying the formula (|T4l) . we obtain the asymptotic solution for the amplitude B 

T3(f\ _ / 0, t < t A ^ B - h, , v 

W \ C B exp [-f ^t + fi^cosfit)] ,t>t^ B + /i, 1 JJ 

where 

(7 = g< ^ A / p 9A->B(tA^B)+fA^B(tA-> B ) (24) 

4iw y (3/2)Fnyi- T 2 /9 

On the second part of the modulation period, the amplitude of the second oscillator is rela- 
tively large, and that of the first oscillator is small. Here we may again regard the coupling as 
unidirectional and include only the effect of the second oscillator on the first one. For the first 
oscillator we write down 



F e 
A + —(7 + sin fit) A = - — C B exp 
2 4uiq 



p 

■— (7* + fi _1 cos fit) 



4tu, 



AeT (25) 



The right-hand part contains two terms. By the linearity of the equation, we will the obtain 
solution as a superposition of two components corresponding to separate contribution of these 
terms: A(t) = A B ^ A {t) + A\(t). 

Let us consider first the solution A B ^ A {t) accounting driving by the second oscillator. Again, 
we may exploit the result for the subsidiary equation setting 

0b-a(*) = (F/2)( 7 + sinfit), f B ^ A (jt) = -(F/2)( 7 t + fi- 1 cosfit), 

K B ^ A = [e/(4u )}C B ^ 

and 

9B-+A® + fB-+A(t) = -FQ- 1 cos fit, 

g' B ^ A (t) + f B _ A (t) = Fsmnt, (27) 
9' B ^ A (t) + f^ A (t)=FQcosQt. 

From the condition of maximum for the function g B ^ A + f B ^ A 

sin flt B ^ A = 0, cosQt B ^ A = — 1 (28) 

we determine a time instant t B ^ A , the neighborhood of which is responsible for excitation of the 
first oscillator. Using formula (|14l) and the relation 

g' B _ A (t B ^ A ) + f^ A (t B ^ A ) = -f3 B ^ A = -Ffi (29) 

we obtain 

J 0, t < t B ^ A - h, 

AB ^ A{t) = { &C B y/%e?/ Q e*j> [-f (t* - IT 1 cos fit)] , t > t B ^ A + h. (30) 

Let us turn now to a solution A\ associated with the second term in the right-hand part of 
the equation (1251) . Now we set 

g' x (t) = (F/2)( 7 + sinfit), / A (t) = 0, K x = -[1 / (4u )}Xe^ (31) 

and 

9x(t) + fx(t) = (l/2)F 7 t - (l/2)Ffi" 1 cos fit, 

g'x(t) + f'x(t) = (l/2)F 7 +(l/2)Fsinfit, (32) 
= (l/2)Ffi cos fit 



7 



Then, we determine t\ corresponding to the maximum of the function g\ + f\ 

smVtt\ = —7, cos fit a = — a/1 — 7 2 . (33) 



Then 

and, from ffT4]) . 



9'{(tx) + fx(tx) = -fix = -(l/2)FQy/T^f (34) 



0, t < t x - h, 

Composing the sum of two solutions (1301) and (1351) . for t > max^^A, £a) we obtain 



AW = <| Ae^ / 4. e3A fa) exp [_£ (7t _ ^-1 CQsQt) j ; t > t A + fe- (35) 



= i^e^ - \e*<> \ exp 



4u v ^ 1 



ry2 



- — (7*- fi" 1 cos fit) 



(36) 



Finally, we evaluate the amplitude of the first oscillator at the end of the period: A T 
A(t n+1 ) = A{t n + T) = A{t n + 2ir/ SI). With the relations (PH), (HMD and ([33ESD this may be 
expressed, via the initial amplitude A n = A(t n ), as 



An+1 ~ 2uJoVFn^~' 



F 

■e 

2 



g (7 arctan -j2=^. + - 7 2 - = J 



^ 2 e 2 v / 2^ r v / l-7 2 ^^ 7 arctan -^j= ~37 arctan +^9-7 2 -3 v /l-7 2 +2-7r7j 
"^i^ov 7 ™ a/9-7 2 

By variable and parameter changes 



(37) 



Ae^ 



_A iV27T£ 2 & (7 arctan -^^arctan -^=+^^-^^+2-2^) 



and 



c = ^ tV2e 2 [7T/(Fn)f/ 2 ^( 7 a r ctan^ + 7arctan^X_ + v ^ 2 - +v ^ 2 - +2 -3 7r7 ) 

16u; 3 </9^yr^ 

the map (137]) is reduced to the canonical form of the complex quadratic map (JT]). 

In Fig. [2] we present diagrams obtained from computations for the map (13T|) and for the 
shortened amplitude equations (J7J). See panels (a) and (b), respectively. Depicted is the plane of 
the complex parameter Xe ltp ; other parameters are N = 10, F = 7, 7 = 0.5, e — 1. Gray color 
corresponds to observation of dynamics with bounded amplitudes, and white to observed escape 
to infinity. The object on the diagram (b) for the set of differential equations (J7|) is evidently 
similar to the Mandelbrot cactus for the complex quadratic map (137|) . Marks 1, 2, 3,... on the 
diagrams indicate leaves, where bounded dynamics of periods T, 2T, 3T, ... take place. Obviously, 
a type of regime is determined by the "germ" signal, which acts in the initial part of the active 
stage of the first oscillator, formed as a composition of the signal from the partner oscillator and 
of the external force. Essential are the phase relations of these two signals, which determine the 
subtle structure of leaves of the "cactus" in the parameter plane. 

Fig. |2Jc), a the point in the parameter plane A = 1.5i, we present a comparison of the time 
dependencies of the amplitudes of the two oscillators obtained by numerical solution of the equa- 
tion (J7J) (gray profiles) and those corresponding to the approximate analytic relations ( TT5|) . (|23|) . 
and ( 136]) (dotted lines). Observe that the numerical and approximate analytic solutions manifest 
good agreement. 
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Figure 2: Mandelbrot set for the map fl37j) (a) and domains of observable bounded periodic dynam- 
ics (gray) in a system of shortened amplitude equations © (b) at parameter values F = 7, 7 = 0.5, 
ujq = 2tt, N — 10, e — 1. The amplitudes of the two oscillators versus time are shown in panel (c) 
for A = 1.51 Gray profiles correspond to numerical solution of the equations ([7]), and dotted 
lines 1-3 designate the approximate analytic solution in accordance with (Tl5|) (1), (123]) (2), ( 136|) (3). 

4 Numerical studies of the basic model system of coupled 
non-autonomous oscillators 

The reduction of the dynamics to the complex quadratic map in the previous section was based 
on a use of the slow-amplitude method (the large- N asymptotics) and of some additional approx- 
imations in the course of derivation of the analytic form of the Poincare map. 

We now return to the original system of coupled non-autonomous oscillators governed by the 
real- value equations ([2]). One may expect that, at least in a coarse sense, the structure of the 
objects in the parameter plane and in the phase space will be similar to that of the Mandelbrot 
and Julia sets of the complex quadratic map. How well the subtler details of the structures are 
present in the original system and in the reduced model, is an interesting and important question. 
In this section, we turn to results of numerical studies of the original system of coupled oscillators. 
In addition, these results may be related to the equivalent set of equations rewritten in terms of 
the complex amplitudes without neglecting the oscillating terms. See (EJ). 

In Fig.|3]we depict the Mandelbrot-like sets obtained from computations for the basic model fl2]). 
The gray areas in the parameter plane corresponding to observation of periodic dynamics in a 
bounded domain of the dynamical variables. The most notable effect in comparison with the slow- 
amplitude approximation consists in a counter-clock-wise rotation of the pictures with a decrease 
of N. Nevertheless, the visible structure of mutual disposition of "leaves" of the "cactuses" 
persists and correspond to that intrinsic to the complex quadratic map (cf. Fig. [1]). Hereafter, to 
avoid redundancy, we will speak of these objects as Mandelbrot cactuses for the coupled oscillator 
system 

For a detailed analysis, let us turn to the case of relatively small N to observe notable deviations 
from the slow-amplitude asymptotic with a possibility of the resolution of the dissimilarity between 
the basic model and the complex analytic map. 

In Fig. S] (panels (a) and (b)) we present pictures of the Mandelbrot cactus obtained in com- 

1 We avoid the term 'Mandelbrot set' because the actual small-scale structure of these objects may be (and is, 
naturally, see below) distinct from that of the classical Mandelbrot set itself. 
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-0.3 A,COS<p 0.3 -0.7 A,COSCp 0.5 -1.25 X COS Cp 0.75 



Figure 3: Mandelbrot cactuses for the system of coupled non-autonomous oscillators (J2J) (or for 
equivalent amplitude equations (jHJ)). Parameter values: (a) N = 80, F = 0.875, e = 0.125, 
(b) N = 40, F = 1.75, e = 0.25, and (c) iV = 20, F = 3.5, e = 0.5. For all cases, co = 2vr 
and 7 = 0.5. Compare with Fig. [2(b) corresponding to the large-iV asymptotic and observe the 
counter-clock-wise rotation of the cactuses under decrease of N. 



putations for the system of coupled non-autonomous oscillators at iV = 10 and other parameters 
uo = 2n, F = 7, 7 = 0.5, e — 1. Except for the approximate right angle turn, the object looks 
similar to those for the complex map and for the shortened amplitude equations. See Fig. [2(a), (b). 
For parameter values marked with dots on the picture of the cactus, we depict diagrams in the 
plane of variables of the first oscillator (x,x/uq), which are analogous to portraits of Julia sets for 
the complex quadratic map in Fig. [TJ To be accurate, we have to note that in the system of coupled 
oscillators we deal with a four-dimensional phase space (x,x/uj ,y,y/2uJo). In the stroboscopic 
Poincare map (as it defined in Section 2), the attractor is placed close to the coordinate plane 
(x,x/uq), but not precisely in it. The basin of attraction is naturally a four-dimensional object. 
The diagrams in panels (g)-(j) of Fig. H] correspond to cross-sections of those four- dimensional 
basins by the plane y = 0, y = for the attractors belonging to a bounded region of the phase 
space. Portraits of the respective attractors in projection onto the plane are shown on panels 
(c)-(f). Dots inside the basins correspond to stroboscopic cross-sections for the attractive peri- 
odic orbits. Fig. [5] illustrates the dynamics of the system on an attractive orbit of period 3T for 
the parameter value Xe itp = —0.2 + 1.5i (it is located inside the "leaf" of the cactus shown with 
magnification in Fig. [1(b)- 

As seen from our work, the model system indeed manifests phenomena known for the complex 
analytic maps, like Mandelbrot and Julia sets, at least on the visible coarse scale level. Moreover, 
it takes place in a wider parameter range than one could expect from the point of view of applica- 
bility of the shortened amplitude equations. Does this similarity extend to the small-scale fractal 
structure of the sets? It appears that this is not the case. The reason for this is the violation of 
complex analyticity as we now illustrate. 

Given an iterative complex analytic map z n+ \ = f(z n ), z = X + iY, one can separate real and 
imaginary parts in the equation and arrive at equivalent description of the dynamics by a real 
two-dimensional map 

X n+1 = U(X n ,Y n ), Y n+1 = V(X n ,Y n ), (40) 

where f(z) = U(X,Y) + iV(X,Y). This is a map of a special kind because the functions must 
satisfy to the Cauchy-Riemann equations 

dU(X,Y) dV(X,Y) dV(X,Y) dU(X,Y) 

OX OY ' OX dY ' 1 ' 

which imply vanishing of a derivative of the function / over the complex conjugate variable: 

df fdU dV\ fdU dV\ 

d?={dx-w) + '{w + ax) =0 - (42) 
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Figure 4: The Mandelbrot cactus (a) and its magnified fragment (b) for the system of coupled non- 
autonomous oscillators (j2J) at uq = 2tt, F = 7, 7 = 0.5, iV = 10, e = 1 in the plane of the complex 
parameter Xe l(p . Gray color corresponds to the presence of attractive periodic orbits. The num- 
bers 1, 2, 3,... designate the period in units of the basic modulation period T = 2n/Q. White color 
corresponds to the absence of attractive periodic motions and escape to infinity. For the parameter 
values marked with dots (panel (a)), portraits of the attractive periodic orbits are shown in projec- 
tion on the plane (x, x/2n) (panels (c) and (f )) together with cross-sections of the basins of attrac- 
tion for these orbits with this coordinate plane ((g) and (j)). Dots inside the basins correspond to 
the stroboscopic cross-sections for those periodic orbits. The panels correspond to the parameter 
sets: A cosy? = 0.5, Asiny? = —0.2 ((c) and (g)), A cosy? = —0.7, Asiny? = 0.4 ((d) and (h)), 
A cos f = —1.7, A sin ip = 0.1 ((e) and (I)), A cosy? = —0.2, A sin tp = 1.5 ((f) and (j)). 
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20.0 




-20.0 



Figure 5: Dynamical variables versus time for the system of coupled non-autonomous oscillators d2J) 
at uq = 2tt, F = 7, 7 = 0.5, e — 1, N — 10, Ae ?v = —0.2 + 1.5i that corresponds to presence of an 
attractive cycle of period 3 in the Poincare map. (The pattern shown is repeated again and again 
with period 3T.) 



Even a small smooth variation of the functions U and V violating the condition of complex 
analyticity implies generically radical changes in the dynamics, e.g. the destruction of the small- 
scale fractal structure of the Mandelbrot set and of intrinsic universal scaling regularities [T7J [T5| 
[IHl [20]. 

Let us look at the complex amplitude equations ([6]), which are equivalent to the original 
equations formulated in real variables. Observe that they contain some terms proportional to 
complex conjugate amplitudes A* and B* . These terms are oscillating, and they disappear under 
the averaging procedure leading to the shortened amplitude equations (j7|). Nevertheless, the 
Poincare map constructed from the exact equation ([6]) inevitably be dependant on the conjugate 
variables, thereby violating the Cauchy-Riemann conditions. This violation may be as small as 
desired in the asymptotic of large N. At relatively small N, the effect of destruction of the 
small-scale structure of the Mandelbrot set becomes observable in computations. 

One tool for an analysis of the degree of violation of the complex analyticity is the com- 
putation of the spectrum of Lyapunov exponents. In particular, for a two-dimensional map 
X n+ i = U(X n ,Y n ), Y n+ i = V(X n ,Y n ) they may be determined via the eigenvalues of the ma- 
trix 

b = a(Xo,ro)a + (Xo,r )a(X 1 ,F 1 )a + (X 1 ,r 1 )...a(X M _ 1 ,r M _ 1 )a + (X M _ 1 ,r M _ 1 ), (43) 

where 

f dU(X,Y)/dX dU(X,Y)/dY\ 

\dV(X,Y)/dX dV(X,Y)/dY)> 1 ] 

In the case of a two-dimensional real map equivalent to an analytic map of one complex variable, 
two Lyapunov exponents must be equal. It may be shown from the Cauchy-Riemann conditions 
that two eigenvalues coincide at any values of parameters and variables. The same is true for the 
Lyapunov exponents expressed as Ai j2 — logAi i 2/2M. 

The non-autonomous system (j2J) possesses four Lyapunov exponents (we exclude the pertur- 
bations associated with a shift along the trajectory, or with a phase of the external driving). To 
compute the Lyapunov exponents we used the Benettin algorithm [2T]. The procedure consists 
in simultaneous numerical solution of the equations ([2]) and a collection of four exemplars of the 
linearized equations for small perturbations: 

x + UqX + F ■ (7 + sin Qt)x = ey sin ui^t, . , 

y + {2uj fy + F ■ (7 - sin Qt)y = 2exx. 

At each period T = 2ir/Q we perform Gram-Schmidt orthogonalization and normalization 
for a set of four vectors x J = {5P, x^/cjq, yi , y^/2oj }, j = 1,...,4. The Lyapunov exponents are 
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estimated as mean rates of growth or decrease of logarithms of the norms of these four vectors: 



1 M 

A i=M^E ln ra, J = l,-,4, (46) 
i=i 

where the norms are evaluated after the orthogonalization but before the normalization. 

Computations show that, depending on the regime, two larger exponents may be negative 
(periodic attractive orbits), positive (chaotic motions) and zero (a border of chaos and quasiperi- 
odic regimes). The other two exponents are always negative in the whole domain of existence of 
bounded dynamical states (i.e. on the Mandelbrot set). In the left column of Fig. E]we present 
charts of the largest Lyapunov exponent on the plane (A sin ip, A cos tp) for three values of N. Gray 
tones from light to dark correspond to variation of the Lyapunov exponent from to — oo. The 
diagram on panel (a) corresponds to iV = 10. Observe that at central parts of the leaves the 
largest Lyapunov exponent becomes large negative, which corresponds to periodic motions of high 
stability. At edges of the leaves, a thin strip of appearance of positive Lyapunov exponent takes 
place (chaos). The picture is similar to that for the map (JT]); see e.g. Ref. [17J. With decrease 
of N, Fig. Etc), distortion of the configuration develops. The leaves lose their round form and 
a wider strip of chaos and quasiperiodicity appears. In the right column of Fig. [61 we depict 
respective charts for the difference of the two larger Lyapunov exponents. In the top diagram this 
difference does not exceed 10~ 2 , which is comparable with numerical errors. However, at smaller 
N (panels (e), (f)) regions of large difference of the exponent appear (black color). It reveals the 
essential deviation from complex analytic dynamics. 

The equality of the pair of larger Lyapunov exponents should be regarded as an indirect 
symptom of complex analytic dynamics. Let us turn to computations aimed at a straightforward 
verification of the Cauchy-Riemann equations (j4ip . For the four-dimensional Poincare map 

x n+1 = F(x n ), (47) 

where x is a vector {x, u,y,v} = {x,x/u ,y,y/2oj }, we produce the following procedure. In the 
course of dynamics of the system on several periods of modulation, we perform numerical solution 
of the equations fl2]) and of the set of equations for four perturbation vectors fj4"5j) with redefinition 
of them at the beginning of each period in accordance with 



^(nT + O) = {x, 0,0,0}, 
X 2 (nT + 0) = {0,u,0,0}, 
£ 3 (nT + 0) = {0,0,^,0}, 
5 4 (nT + 0) = {0,0, 0,v}. 

At the end of each period, we compose a matrix 

/5 n (nT-0)/x 5t 21 (nT-0)/u 5i 31 (nT-0)/y x 41 (riT - 0)/v \ 

x 12 (nT-0)/x 5 22 (nT-0)/w X 32 (nT-0)/y x. 42 (nT-0)/v 

5 13 (nT-0)/x 5 23 (nT-0)/w 5i 33 (nT-0)/y 5c 43 (nT-0)/v 

y x 14 (nT - 0) jx 5 24 (nT - 0)/u x 34 (nT - 0) /y 5 44 (nT -0)/v J 



(48) 



A 



(49) 



In the case of exact fulfillment of the conditions of analyticity of the map (T4T|) the elements must 
satisfy 

A 11 = A 22 A 12 = -A 21 A 31 = A 42 , A 32 = -A 41 
A 13 = A 24 A 14 = -A 23 A 33 = A 



13 A - I AN a 2:? a \ ! i a 13 a I (50) 



Alternatively, it is convenient to consider derivatives of the complex functions F\ and F2, defined 
as components of the vector function ( 14T|) over the conjugate variables p* = x — iu and q* = y — iv, 
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Figure 6: Charts of the largest Lyapunov exponent (a-c) and for the difference of the two larger 
exponents (d-f) for the system of coupled non-autonomous oscillators ^2} at three values of the 
ratio of periods of modulation and of basic oscillators: N = 10 ((a), (d)), N = 6 ((b), (e)), 
iV = 3 ((c), (f)). Other parameters are as in Fig. |H Uniform gray color means area of unstable 
dynamics (typically divergence to infinity). The legend for gray scales is shown at the bottom of 
each column. Black color on the diagrams in the left column also corresponds to chaotic dynamics 
in a bounded domain. 
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Figure 7: The logarithm of the absolute value of the derivative of one of the functions determining 
the two-dimensional complex Poincare map p n +\ = Fi(p n ,q n ), q n+ i = -F 2 (p n ,g n ) with respect to 
the complex conjugate variable. The complex variables are expressed via real physical variables 
of the system (j2J) as p = x + ix/uo, q = y + iy/2ujQ. The derivative was evaluated at the origin 
x = 0, x = 0, y = 0, y = for A = 0. Other parameters are as in Fig. HI 



from which one can diagnose the presence or absence of the non-analyticity: 

fi = (A 11 - A 22 ) + *(A 12 + A 21 ), ff = (A 31 
§| = (A 13 -A 24 ) + *(A 14 + A 23 ), || = (A 33 



A 42 ) + z(A 32 + A 41 ), 
A 44 ) + i(A 43 + A 34 ). 



(51) 



In Fig. we present a plot of the logarithm of the absolute value of the derivative \dF\/dp*\ 
at the origin (x(0) = (0,0,0,0)) with A = versus the parameter N, that is the dimensionless 
period of modulation. This is a quantifier or a degree of non-analyticity of the function under 
consideration. The data are approximated by a straight line with slope —1.86. It means that the 
degree of non-analyticity determined at the origin manifests exponential decay with growth of N 
as Fi(p,p*,q,q*) ~ e~ im p* . 

Figure shows diagrams illustrating the distribution of the maximal values of ratios 



8F-, 



dp* 

dF 2 
dp* 



dF, 



dp 

dF 2 
dp 



do 



8Fi 



dq* 

dF 2 
dq* 



dF, 



dq 

dF 2 
dq 



(52) 



on the parameter plane \e tip . The values are determined along the orbits starting from the origin of 
length equal to 100 modulation periods. Gray scales from light to dark correspond to growth of d, 
i.e. to increase of degree of deflection of the dynamics from the pure complex analytical. Observe 
that at small iV (panel (c)), wide regions of violation of the Cauchy-Riemann conditions take place, 
where the ratios of derivatives in respect to complex variable and the conjugate variable exceed 1. 
For larger N, when the cactuses look similar to the Mandelbrot set (panels (a) and (b)), analogous 
violations occur only on the small-scale details of the Mandelbrot-like structure (see Fig. 

It is interesting to evaluate degree of violation of the analyticity globally in the phase space of 
the system or, at least, in domains including the basins of attraction. In Fig. [10] we present data 
of computations aimed at estimating the maximal absolute value for the derivative ratios fl52l) in 
the plane (x,u). Observe that with a decrease of the parameter N the picture becomes darker. 
(The color coding is assumed the same as in Fig.[S](a)-(c).) In Fig.[TT]we present plots for maximal 
values of di (i = 1, ...,4) determined in a domain of a four- dimensional cube in the phase space 
containing the attraction basins. (Computations were produced at nodes of four-dimensional grid 
of size 50 x 50 x 50 x 50.) Observe that at small N the derivatives over the conjugate variables 
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Figure 8: The distributions of the maximal (a-c) and average (d-f) ratios of derivatives f)49p over 
the plane of the complex parameter Xe lip . The values are determined along the orbits starting from 
the origin, of length equal to 100 modulation periods. Gray scales from light to dark correspond 
to increase of degree of deviation of the dynamics from the pure complex analytical. The diagrams 
are drawn for iV = 10 (a,d), N = 6 (b,e), N = 3 (c, /). Other parameters are as in Fig. |H 
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Figure 10: The distributions of the maximal ratios of derivatives (1491) over the plane of the complex 
parameter Xe l(p determined for one modulation period. The diagrams are drawn for N = 3 (a), 
N = 6 (b), N = 10 (c), N = 15 (d). Other parameters are as in Fig. |H 




Figure 11: The logarithm of the maximal (left panel) and average (right panel) value of the ratios 
of derivatives d computed on one iteration of the Poincare map over nodes of an array of size 
50 x 50 x 50 x 50 on a domain of a four-dimensional cube in the phase space containing the basins 
of attraction. The diagram is drawn for A = 0, uq = 2tt, F = 7, 7 = 0.5, e — 1. 



become rather large in comparison with derivatives in respect to the main complex variable; for 
N < 4 they may be larger even by many times. 

As follows from our computations, the system of coupled non-autonomous oscillators indeed 
demonstrates dynamics roughly corresponding to that in the complex analytic map. The degree 
of the correspondence is determined by parameter N, that is a ratio of modulation period to 
period of basic oscillations. With decrease of N, the type of dynamics is changed gradually; the 
destruction of the picture associated with the complex analytic dynamics starts from small-scale 
details of the visible structure of the Mandelbrot set. 



5 Complex analytic dynamics in coupled self-sustained os- 
cillators 

Let us modify original system (j2J) in following way 

x + UqX + [F ■ (7 + sin fit) + 5x 2 ] x = eysinuot + Xsm(uot + <p), , . 

y + (2co ) 2 y + [F ■ (7 - sin fit) + 5y 2 ]y = ex 2 , ^ } 
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Figure 12: Senior Lyapunov exponent chart, demonstrating Mandelbrot set at the parameter 
plane of the system of coupled self- sustained oscillators (|5"3"|) at 5 = 0.00001, F = 2.0, 7 = 0.5, 
e = 1.5, N = 32, uj — 2it (a). Portraits of the associated with Mandelbrot set attractor of period 2 
(b), corresponded plot of x and y dynamical variables versus time (d) and multistable alternative 
attractor (c) at Acos<£> = —0.03, Asin<£> = 0.045. Homogenous gray color on panel (a) corresponds 
to an escape of the trajectories far away from origin and convergence to alternative attractor. 

This is a system of two coupled non-autonomous alternately excited self- sustained van der Pol 
oscillators. With 5 = it comes to investigated system (j2J). With A = it comes to the proposed 
in [10] model, manifesting hyperbolic chaos and Smale- Williams attractor. 

Fig. 12 (a) demonstrates the chart of the parameter plane (A cos ip , Asin^), at which the struc- 
ture similar to the Mandelbrot set is visible. The regions around the Mandelbrot set, corresponds 
to escaping of the trajectories far away from the vicinity of origin - domain at which attractors, 
associated with Mandelbrot set can exist. However, in contrast with the system (J2J), where tra- 
jectories diverge and goes to infinity, in coupled self-sustained oscillators due to nonlinear terms, 
proportional to xx 2 and yy 2 , trajectories comes to mulltistable periodic or chaotic attractor. Co- 
existence of two attractors - associated with Mandelbrot set, and alternative one is shown in 
Fig. 12 (b) and 12 (c) respectively. The factor 5 of nonlinear terms xx 2 and yy 2 drives the charac- 
teristic size of alternative attractor. It shrinks with increasing S and can destroy "complex analytic 
dynamics" near origin. 
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6 Conclusion 



In this paper, a system was proposed that consists of two coupled alternately excited oscillators 
with a turn-by-turn transfer of the excitation from one to another, accompanied with appropriate 
nonlinear transformation of the complex amplitude of the oscillations in the course of the process. 
This system obviously allows realization as a physical object, e.g. as an electronic device analogous 
to that described in Ref. [TD]. Analytic consideration shows that the Poincare map for the system 
corresponds to a definite approximation to a complex quadratic map. Numerical studies confirm 
the presence of the expected phenomena intrinsic to iterative complex analytic maps, namely 
the Mandelbrot set and Julia sets, at least up to a definite level of resolution of the fractal-like 
structures. Analysis of the violation of the applicability of the approximation corresponding to 
the complex quadratic map revealed several effects. One is the rotation of the Mandelbrot- like set 
in the complex parameter plane. Another is the destruction of the small-scale fractal structure 
under a decrease of the parameter representing the ratio of the modulation period to the period 
of basic oscillations. 
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